clc; clear; close all;

% 参数设置
mu = 3.5;           
z0 = -1;            
k_values = [0.15, 0.2, 0.25, 0.3];  
nx = 100;           % 定义清晰度           
ny = nx;           
x_range = linspace(0, 1, nx);  
y_range = linspace(0, 1, ny);  

max_iter = 1000;    % 总迭代步数
transient = 500;    % 瞬态步数
tol = 1e-3;         % 周期检测容差
max_range = 10;     % 系统最大范围，超过此范围认为发散

% 定义循环缓冲区大小（用户可自定义，需大于最大周期检测 p 的最大值）
buff_size = 21;  % 例如：21 表示可检测周期 1 到 20

% 用于保存分类结果的矩阵，存储每个 (x0,y0) 点的状态代码
% 状态代码定义：
% 1  : SP, 2 : P2, 4 : P4, 8 : P8, 9 : MP, 12 : QP, 13 : CH, 14 : HC, 15 : UB
basin_results = zeros(ny, nx, length(k_values));

% 定义自定义 colormap（多段渐变）
cmap_base = [...
    75,   0, 130;   % Indigo (SP)
    0, 255, 0;      % Green (P2)
    255, 215, 0;    % Gold (P4)
    255, 140, 0;    % DarkOrange (P8)
    255, 0, 0;      % Red (MP)
    0, 255, 255;    % Cyan (QP)
    255, 0, 255;    % Magenta (CH)
    255, 165, 0;    % Orange (HC)
    128, 128, 128;  % Gray (UB)
];

cmap_base = cmap_base / 255;  % Normalize
cmap_custom = interp1(linspace(0,1,size(cmap_base,1)), cmap_base, linspace(0,1,256));

for idxK = 1:length(k_values)
    k = k_values(idxK);
    fprintf('Processing k = %g\n', k);
    
    % 遍历 (x0,y0) 网格
    for ix = 1:nx
        for iy = 1:ny
            % 初始化当前点的状态
            state = [x_range(ix), y_range(iy), z0];
            
            % 进行 transient 步（不存储轨迹）
            for n = 1:transient
                [dx, dy, dz] = mclm(state, mu, k);
                state = [dx, dy, dz];
            end
            
            % 利用循环缓冲区保存最后 buff_size 次迭代的 (x,y) 状态
            buffer = zeros(buff_size, 2);
            total_eval = max_iter - transient;  % 正式迭代步数
            for n = 1:total_eval
                [dx, dy, dz] = mclm(state, mu, k);
                state = [dx, dy, dz];
                idx_buf = mod(n-1, buff_size) + 1;  % 当前缓冲区索引
                buffer(idx_buf, :) = state(1:2);
            end
            
            % 计算缓冲区中最后一个状态的索引
            last_index = mod(total_eval-1, buff_size) + 1;
            % 调用分类函数对当前点进行状态分类
            code = classify_5(buffer, last_index, buff_size, tol, max_range, mu, k);
            basin_results(iy, ix, idxK) = code;
        end
    end
    
    % 用 subplot 绘制当前 k 值下的吸引盆图
    subplot(2,2,idxK);
    imagesc(x_range, y_range, basin_results(:,:,idxK));
    set(gca, 'YDir', 'normal');
    colormap(gca, cmap_custom);
    colorbar;
    title(['k = ', num2str(k), ', \mu = ', num2str(mu), ', z_0 = ', num2str(z0)]);
    xlabel('x_0');
    ylabel('y_0');
end